# Replication of all figures in the paper and most figures in the appendix
# These are partial replication files, as some of the data is pre-processed
# Estimates of the position of presidents and parties are available as a 
# data file in .csv and .RData formats
# Please contact cesar.zucco@fgv.br for additional replication information

# (1) Place all datafiles and code in the same folder
# data_for_polarization.RData
# bls9_estimates_medianlegislator.RData
# bls9_estimates_partiespresidents_long.RData

# (2) Obtain the PREPPS data from  and place it in the same folder
# PREPPS data can be obtained from http://wiesehomeier.net/prepps/

# (3) Code will generate the following output (in the same folder)
# Figure 1
# Figure 2 
# Figure 3 
# Figure D1
# Figure E1

library(plyr)
library(tidyverse)

## Figure 1 ######################################################

# Declare the function that produces the plot
plot.party.year<-function(x=to.plot,pe.string="ideo",
                          se.string="ideo.se",
                          draw.arrows=FALSE,
                          arrow.col=2,label.years=FALSE){
  require(plyr)
  require(reshape)
  #x are point estimates, y are dispersion/uncertainty estimates
  par(mfrow=c(1,1), mar=c(2,6.5,1,1.5))
  pe <- x[,pe.string]
  se <- x[,se.string]
  #use last position of parties to order parties
  tmp <- ddply(x,"party.or.pres",function(y) y[nrow(y),pe.string])
  parties <-  tmp$party.or.pres[sort(tmp$V1,index.return=T)$ix]
  #basic plot
  plot(range(na.omit(c(pe-1.64*se,pe+1.64*se))),
       c(1,length(parties)), type="n",
       yaxt="n",xaxt="n",ylab="",
       xlab="",bty="n")
  axis(2,at=1:length(parties),
       labels=parties,srt=90,las=2,tick=FALSE,
       cex.axis=1.5,line=-.3)
  mtext(expression("" %<-% Left),side=1,line=.5, adj = 0,cex=1.5)
  mtext(expression(Right %->% ""),side=1,line=.5, adj = 1,cex=1.5)
  for(p in 1:length(parties)){
    y2<-p+0.5
    y1<-p-0.5
    polygon(x=c(-10,10,10,-10),y=c(y1,y1,y2,y2),
            col=if(round(p/2)==p/2){gray(0.8)}else{gray(0.6)},border=NA)
    by.party <- subset(x,party.or.pres==parties[p])
    party.years <- sort(unique(x$year))
    struc <- data.frame(year=party.years, 
                        ys = seq(p+.25,p-0.25,length=length(unique(x$year))))
    by.party <- merge(by.party,struc,by="year",all.y=T)
    pe <- by.party[,pe.string]
    se <- by.party[,se.string]
    ys <- by.party$ys
    if(label.years==TRUE){
      axis(4,at=ys[c(1,length(ys))],
           labels=party.years[c(1,length(ys))],srt=90,
           las=2,tick=FALSE,
           cex.axis=0.7,line=-.9)
      abline(h=ys,col=gray(0.9),lwd=0.3)
    }
    points(pe,ys,pch=20,cex=1.3)
    segments(x0=pe-1.64*se,
             x1=pe+1.64*se,
             y0=ys,y1=ys,lwd=2)
    if(draw.arrows==TRUE){
      arrows(x0=pe[1],x1=pe[4],
             y0=ys[1],y1=ys[4],
             col=arrow.col, 
             code=2,lwd=2,length=0.05,cex=0.5)
      arrows(x0=pe[4],x1=pe[6],
             y0=ys[4],y1=ys[6],
             col=arrow.col,#gray(0.6),
             code=2,lwd=2,length=0.05,cex=0.5)
      arrows(x0=pe[6],x1=pe[8],
             y0=ys[6],y1=ys[8],
             col=arrow.col, 
             code=2,lwd=2,length=0.05,cex=0.5)
      arrows(x0=pe[8],x1=pe[length(pe)],
             y0=ys[8],y1=ys[length(ys)],
             col=arrow.col, 
             code=2,lwd=2,length=0.1,cex=0.5) 	
    }#end arrows
  }
}

#Load output from rescaling procedure
load('~/Dropbox/Data/BLS9-Rescale/bls9_estimates_partiespresidents_long.RData')

#Select parties to report
to.plot <- subset(long.table,
                  is.element(party.or.pres,c("CID","PCDOB","PP","PDT",
                                     "DEM","PL","MDB","PSB","PSDB",
                                     "PT","PTB","PSOL","PSD","REP",
                                     "NOVO","PSL")))

#Generate jpeg plot in grayscale 
jpeg("fig_1partiesselected_arrow.jpg"
    ,width=8,height=8,units="in",res=300)
plot.party.year(to.plot,label.years=T,draw.arrows=T,arrow.col=gray(0.4))
dev.off()



## Figure 2 ######################################################

# Declare function that produces the plot
vertical.plot<-function(x=to.plot,
                        pe.string="ideo",
                        se.string="ideo.se",
                        parties=NULL,
                        medians=medw){#
  #x are point estimates, y are dispersion estimates, m are legislator median est
  all.years <- sort(unique(x$year))
  if(is.null(parties)){parties<-levels(factor(x$party.or.pres))}
  pe <- x[,pe.string]
  se <- x[,se.string]
  par(mfrow=c(1,1))
  plot(medians,length(all.years):1
       ,xlim=range(c(pe+1.64*se,
                     pe-1.64*se))
       ,xaxt="n",lty=2,yaxt="n",xlab="",ylab="",bty="n")
  mtext(expression("" %<-% Left),side=1,line=.8, adj = 0,cex=1.5)
  mtext(expression(Right %->% ""),side=1,line=1.1, adj = 1,cex=1.5)
  the.labels <- paste(all.years,"\nBLS ",1:length(all.years),sep="")
  axis(2,at=length(all.years):1, labels=the.labels
       , cex.axis=1.3, las=1)
  abline(h=1:length(all.years),lty=3,col=gray(0.7))
  lines(medians,length(all.years):1,col=1,lty=2,lwd=2)
  cols <- gray.colors(length(parties),end = 0.7)
  for(p in 1:length(parties)){
    by.party <- subset(x,party.or.pres==parties[p])
    years.party <- sort(unique(by.party$year))
    the.years <- which(is.element(all.years,years.party))
    pe <- by.party[,pe.string]
    se <- by.party[,se.string]
    points(pe,rev(the.years),pch=20,cex=2,col=cols[p])
    lines(pe,rev(the.years),lty=2,col=cols[p])
    for(i in nrow(by.party):1){
      lines(c(pe[i]-1.64*se[i],pe[i]+1.64*se[i]), 
            c(rev(the.years)[i],rev(the.years)[i]),
            lwd=2,col=cols[p])  }
  }   
  tmp <- subset(to.plot,is.element(to.plot$party.or.pres,parties)&to.plot$year==1990)
  axis(side=3,at=tmp[,pe.string],labels=tmp[,"party.or.pres"],tick=F,line=-1.5,cex.axis=.9)
  axis(side=1,at=medians[length(medians)],labels="Median",tick=F,line=-1.5)
}

#Load output from rescaling procedure
load("bls9_estimates_partiespresidents_long.RData")
to.plot <- long.table

#Load estimates for "median" legislator
load("bls9_estimates_medianlegislator.RData")

jpeg("fig_2presidentsandtheirparties.jpg",
      width=8,height=8,units="in",res=300)
#Produce the figure basic figure
par(mar=c(6,4,2,2))
vertical.plot(to.plot,parties=c("")
              ,medians= the.medians$med.legis.w )

#Add parties
the.pt <- which(to.plot$party.or.pres=="PT"&is.element(to.plot$year,c(2005,2009,2013,2017)))
points(to.plot$ideo[the.pt],5:2,pch=24,col=1,bg=c(1))
segments(x0=to.plot$ideo[the.pt]+1.64*to.plot$ideo.se[the.pt],
         x1=to.plot$ideo[the.pt]-1.64*to.plot$ideo.se[the.pt],
         y0=5:2,y1=5:2,col=1,lwd=1)
the.mdb <- which(to.plot$party.or.pres=="MDB"&is.element(to.plot$year,c(1990,1993,2017)))
points(to.plot$ideo[the.mdb],c(9,8,2),pch=24,col=c(1,gray(0.5),gray(0.5)),bg=c(1,gray(0.5),gray(0.5)))
segments(x0=to.plot$ideo[the.mdb]+1.64*to.plot$ideo.se[the.mdb],
         x1=to.plot$ideo[the.mdb]-1.64*to.plot$ideo.se[the.mdb],
         y0=c(9,8,2),y1=c(9,8,2),col=c(1,gray(0.5),gray(0.5)),lwd=1)
the.psdb <- which(to.plot$party.or.pres=="PSDB" & is.element(to.plot$year,c(1997,2001)))
points(to.plot$ideo[the.psdb],7:6,pch=24,col=1,bg=1)
segments(x0=to.plot$ideo[the.psdb]+1.64*to.plot$ideo.se[the.psdb],
         x1=to.plot$ideo[the.psdb]-1.64*to.plot$ideo.se[the.psdb],
         y0=7:6,y1=7:6,col=1,lwd=1)
the.prn <- which(to.plot$party.or.pres=="PRN" & is.element(to.plot$year,c(1993)))
points(to.plot$ideo[the.prn],8,pch=24,col=1,bg=1)
segments(x0=to.plot$ideo[the.prn]+1.64*to.plot$ideo.se[the.prn],
         x1=to.plot$ideo[the.prn]-1.64*to.plot$ideo.se[the.prn],
         y0=8,y1=8,col=1,lwd=1)
the.psl <- which(to.plot$party.or.pres=="PSL"&is.element(to.plot$year,c(2021)))
points(to.plot$ideo[the.psl],1,pch=24,col=1,bg=1)
segments(x0=to.plot$ideo[the.psl]+1.64*to.plot$ideo.se[the.psl],
         x1=to.plot$ideo[the.psl]-1.64*to.plot$ideo.se[the.psl],
         y0=1,y1=1,col=1,lwd=1)


#Add presidents
the.sarney <- which(to.plot$party.or.pres=="SARNEY"&to.plot$year==2013)
text(to.plot$ideo[the.sarney],9,pos=1,labels="Sarney")
segments(x0=to.plot$ideo[the.sarney]+1.64*to.plot$ideo.se[the.sarney],
         x1=to.plot$ideo[the.sarney]-1.64*to.plot$ideo.se[the.sarney],
         y0=9,y1=9,col=1,lwd=2)
points(to.plot$ideo[the.sarney],9,pch=22,col=1,bg="white")
the.collor <- which(to.plot$party.or.pres=="COLLOR"&to.plot$year==2013)
text(to.plot$ideo[the.collor],8,pos=1,labels="Collor*")
segments(x0=to.plot$ideo[the.collor]+1.64*to.plot$ideo.se[the.collor],
         x1=to.plot$ideo[the.collor]-1.64*to.plot$ideo.se[the.collor],
         y0=8,y1=8,col=1,lwd=2)
points(to.plot$ideo[the.collor],8,pch=22,col=1,bg="white")
the.itamar <- which(to.plot$party.or.pres=="ITAMAR"&to.plot$year==2013)
text(to.plot$ideo[the.itamar],8,pos=1,labels="Itamar",col=gray(0.5))
segments(x0=to.plot$ideo[the.itamar]+1.64*to.plot$ideo.se[the.itamar],
         x1=to.plot$ideo[the.itamar]-1.64*to.plot$ideo.se[the.itamar],
         y0=8,y1=8,col=gray(0.5),lwd=2)
points(to.plot$ideo[the.itamar],8,pch=22,col=gray(0.5),bg="white")
the.fhc <- which(to.plot$party.or.pres=="FHC"&to.plot$year==2013)
text(to.plot$ideo[the.fhc],7,pos=1,labels="FHC")
segments(x0=to.plot$ideo[the.fhc]+1.64*to.plot$ideo.se[the.fhc],
         x1=to.plot$ideo[the.fhc]-1.64*to.plot$ideo.se[the.fhc],
         y0=7,y1=7,col=1,lwd=2)
points(to.plot$ideo[the.fhc],7,pch=22,col=1,bg="white")
the.fhc2 <- which(to.plot$party.or.pres=="FHC"&to.plot$year==2013)
text(to.plot$ideo[the.fhc2],6,pos=1,labels="FHC")
segments(x0=to.plot$ideo[the.fhc2]+1.64*to.plot$ideo.se[the.fhc2],
         x1=to.plot$ideo[the.fhc2]-1.64*to.plot$ideo.se[the.fhc2],
         y0=6,y1=6,col=1,lwd=2)
points(to.plot$ideo[the.fhc2],6,pch=22,col=1,bg="white")
the.lula <- which(to.plot$party.or.pres=="LULA"&to.plot$year==2013)
text(to.plot$ideo[the.lula],5,pos=1,labels="Lula")
segments(x0=to.plot$ideo[the.lula]+1.64*to.plot$ideo.se[the.lula],
         x1=to.plot$ideo[the.lula]-1.64*to.plot$ideo.se[the.lula],
         y0=5,y1=5,col=1,lwd=2)
points(to.plot$ideo[the.lula],5,pch=22,col=1,bg="white")
the.lula2 <- which(to.plot$party.or.pres=="LULA"&to.plot$year==2013)
text(to.plot$ideo[the.lula2],4,pos=1,labels="Lula")
segments(x0=to.plot$ideo[the.lula2]+1.64*to.plot$ideo.se[the.lula2],
         x1=to.plot$ideo[the.lula2]-1.64*to.plot$ideo.se[the.lula2],
         y0=4,y1=4,col=1,lwd=2)
points(to.plot$ideo[the.lula2],4,pch=22,col=1,bg="white")
the.dilma <- which(to.plot$party.or.pres=="DILMA"&to.plot$year==2013)
points(to.plot$ideo[the.dilma],3,pch=22,col=1,bg=1)
text(to.plot$ideo[the.dilma],3,pos=1,labels="Dilma")
segments(x0=to.plot$ideo[the.dilma]+1.64*to.plot$ideo.se[the.dilma],
         x1=to.plot$ideo[the.dilma]-1.64*to.plot$ideo.se[the.dilma],
         y0=3,y1=3,col=1,lwd=2)
the.dilma2 <- which(to.plot$party.or.pres=="DILMA"&to.plot$year==2017)
points(to.plot$ideo[the.dilma2],2,pch=22,col=1,bg=1)
text(to.plot$ideo[the.dilma2],2,pos=1,labels="Dilma*")
segments(x0=to.plot$ideo[the.dilma2]+1.64*to.plot$ideo.se[the.dilma2],
         x1=to.plot$ideo[the.dilma2]-1.64*to.plot$ideo.se[the.dilma2],
         y0=2,y1=2,col=1,lwd=2)
the.temer <- which(to.plot$party.or.pres=="TEMER"&to.plot$year==2017)
points(to.plot$ideo[the.temer ],2,pch=22,col=gray(0.5),bg=gray(0.5))
text(to.plot$ideo[the.temer ],2,pos=1,labels="Temer",col=gray(0.5))
segments(x0=to.plot$ideo[the.temer ]+1.64*to.plot$ideo.se[the.temer ],
         x1=to.plot$ideo[the.temer ]-1.64*to.plot$ideo.se[the.temer ],
         y0=2,y1=2,col=gray(0.5),lwd=2)
the.bozo <- which(to.plot$party.or.pres=="BOLSONARO"&to.plot$year==2021)
points(to.plot$ideo[the.bozo],1,pch=22,col=1,bg=1)
text(to.plot$ideo[the.bozo],1,pos=1,labels="Bolsonaro",xpd=NA)
segments(x0=to.plot$ideo[the.bozo]+1.64*to.plot$ideo.se[the.bozo],
         x1=to.plot$ideo[the.bozo]-1.64*to.plot$ideo.se[the.bozo],
         y0=1,y1=1,col=1,lwd=2)

legend(x="bottom",pch=c(22,22,24),pt.bg=c("white","black","black")
       ,legend=c("President (retrospective)",
                 "President (contemporaneous)",
                 "President's party")
       ,inset=-0.2
       ,xpd=NA, bty="n")
dev.off()


#### Figure 3 #####

#Load data on parties (assembled from many sources)
load("data_for_polarization.RData")

#Get ideological placement through 2021 and merge in 
load('~/Dropbox/Data/BLS9-Rescale/bls9_estimates_parties_long.RData')
load('~/Dropbox/Data/BLS9-Rescale/bls9_estimates_partiespresidents_long.RData')


id <- party.table[,1:3]
names(id) <- c("wave","party","ideo")
id$year <- ifelse(id$wave==1990,1986,
                  ifelse(id$wave==1993,1990,as.numeric(id$wave)-3))
id$party <- car::recode(id$party,"'PCDOB'='PCdoB'")
id$party <- car::recode(id$party,"'CID'='PPS'")#backwards compatibility with legislative data
id$party <- car::recode(id$party,"'PL'='PR'")
id$party <- car::recode(id$party,"'MDB'='PMDB'")
id$party <- car::recode(id$party,"'NOVO'='Novo'")
id$party <- car::recode(id$party,"'REP'='PRB'")
id$party <- car::recode(id$party,"'REDE'='Rede'")

  # Some guesswork (See paper for details)
  # For 57th legislature: use estimates based on BLS9, but:
  # combining DEM and PSL for UNIAO
  # combining BOLSONARO and PL for PL 
  id2 <- subset(id,wave==2021)
  id2$wave <- id2$year <-2022
  weight.psl <- 0.0926/(0.0926+0.0556) #relative sizes at time of 2022 survey
    weight.dem <- 0.055/(0.0926+0.0556) #relative sizes at time of 2022 survey
  id2$ideo[which(id2$party=="DEM")] <- id2$ideo[which(id2$party=="DEM")] * weight.dem +
  id2$ideo[which(id2$party=="PSL")] * weight.psl
  id2 <- id2[-which(id2$party=="PSL"),] #remove PSL
  weight.pl <- 45 / (112)
  weight.b <-  (112-45) / (112)
  id2$ideo[which(id2$party=="PR")] <- long.table$ideo[which(long.table$party.or.pres=="BOLSONARO")] * weight.b  +
  long.table$ideo[which(long.table$party.or.pres=="PL"&long.table$year=="2021")] * weight.pl
  id <- rbind(id,id2)
  rm(party.table)
  id$party <- as.character(id$party)
  id$party[which(id$party=="DEM"&id$year==2022)] <- "UNIAO"
  id$party <- factor(id$party)
  #end guess work

id$party <- factor(id$party)
ddid <- merge(ddp,id,by=c("party","year"),all.x=T)

# manually enter seat shares in 1984
ddid$ss[which(ddid$year==1982)]<-c(4.8,41.8,49.1,1.7,2.7)/100
ddid$wave[which(ddid$year==1982)]<- 1984

### Compute polarization indices #####

### Declare polarization index function ###
polmaad <- function(x,type="vs"){
  x <- na.omit(x);
  if(type=="vs"){x$s<-x$vs}else{x$s<-x$ss}
  x$s <- x$s/sum(x$s) #compute wights: vs relative to parties for which we have ideo
  insum <- rep(NA,nrow(x))
  for(i in 1:nrow(x)){
    insum[i] <- x$s[i] * sum(x$s[-1] * abs(x$ideo[i]-x$ideo[-i]))
  }
  mad <- sum(insum)
  return(mad)}	

# Stats with SeatShares
polms <-  ddply(ddid, ~ year, polmaad,type="ss")
names(polms)[2]<-"maad"

# Seats share for which we have ideo data
ss <- ddply(ddid, ~ year, function(x){c(tss=100*sum(x$ss[which(is.na(x$ideo)==F)]),
                                        n=length(x$ss[which(is.na(x$ideo)==F)]))})
poltabss <- merge(ss,polms,by="year")

# Get presidential positions to compute distance from presidents to parties
pres <- subset(long.table,is.element(party.or.pres
                                     ,c("SARNEY","COLLOR","ITAMAR","LULA","DILMA","TEMER","BOLSONARO")))

# Create an object with the president and BLS year that she/he was serving
pres <- subset(long.table,
               party.or.pres=="SARNEY"&year==2013|
                 party.or.pres=="COLLOR"&year==2013|
                 party.or.pres=="ITAMAR"&year==2013|
                 party.or.pres=="FHC"&year==2013|
                 party.or.pres=="LULA"&year==2013|
                 party.or.pres=="DILMA"&year==2013|
                 party.or.pres=="DILMA"&year==2017|
                 party.or.pres=="TEMER"&year==2017|
                 party.or.pres=="BOLSONARO"&year==2021)
pres <- rbind(pres,pres[3,],pres[5,]) #repeat two term presidents FHC and Lula
pres$year <- c(1993,2013,1997,1993,2005,1990,2017,2017,2021,2001,2009)
pres <- pres[order(pres$year),]
rownames(pres) <- NULL

## Add Lula 57th legis for "approximate" calculations
lula.57 <- subset(long.table,party.or.pres=="LULA"&year==2017)
lula.57[1] <- 2022
pres <- rbind(pres,lula.57)

# Use BLS years in which presidents were in office
prid <- merge(ddid,pres[,1:4]
              ,by.x="wave",by.y="year"
              ,suffixes=c("",".pres"))
pd.maad <- function(x,type="vs"){
  x <- na.omit(x);
  if(type=="vs"){x$s<-x$vs}else{x$s<-x$ss}
  x$s <- x$s/sum(x$s) #compute vs relative to parties for which we have ideo
  mad <- sum(x$s *abs(x$ideo-x$ideo.pres))
  return(mad)}	

pres.maad <-  ddply(prid, ~ wave + party.or.pres, pd.maad)
names(pres.maad) <- c("wave","pres","maad")

## Size of centrist parties
ddid$center <- ddid$ideo> -0.25 & ddid$ideo < 0.25
the.center <- subset(ddid,center==T)
center <- as.data.frame(as.matrix(by(the.center$ss,the.center$year,sum)))
names(center) <- "center"
poltabss <- merge(poltabss,center,by.x="year",by.y=0)
poltabss$legis <- 48:57

### Plot the figure in the paper ###
jpeg(file="fig_3threewayfigure.jpg",
     width=8,height=6,units="in",res=300)
par(mar=c(6,4,2,4))
tmp <- barplot(by(the.center$ss,the.center$year,sum)
               ,yaxt="n",ylab="",xaxt="n"
               ,density=c(rep(NA,9),20)
               ,xlim=c(0,12))
axis(side=4,at=c(0,0.1,0.2,0.3,0.4))
mtext("Centrist Share of Congress (bars)",side=4,line=2.5)
poltabss$y <- tmp[,1]

par(new=T)
nn <- nrow(poltabss)
plot(poltabss$y,poltabss$maad,type="l"
     ,xaxt="n",xlab="",bty="n",ylim=c(0.4,0.8)
     ,xlim=c(0,nn+2)
     ,ylab="Mean Weighted Distances (line and dots)")
lines(poltabss$y[(nn-1):nn],
      poltabss$pol.maad[(nn-1):nn],
      lty=2,col=gray(0.9))
axis(side=1,at=poltabss$y,labels=paste("Legis\n",poltabss$legis),padj=.5)

#Modify presidential object to include in figure
pres.maad$y <- c(tmp[1,1],tmp[2,1]+c(-.3,.1)
                 ,tmp[3:7,1],tmp[8,1]+c(-.3,.1),tmp[9:10,1])
points(pres.maad$y,pres.maad$maad,pch=21
       ,bg=c(rep("black",length(pres.maad$y)-1),gray(0.9)))
text(pres.maad$y,pres.maad$maad,cex=0.7,labels=pres.maad$pres
     ,pos=3,xpd=NA)

legend(x="bottom"
       ,legend=c("Distance between parties","Distance between president and parties")
       ,lty=c(1,NA)
       ,pch=c(NA,21)
       ,pt.bg=c(NA,"black")
       ,bty="n",horiz=T
       ,inset=-0.25,xpd=NA)
dev.off()


### Plot with multiple definitions of "center" in appendix 

## Several definitions of center
all.centers <- poltabss
for(i in seq(0.20,0.45,by=0.05)){
  ddid$center <- ddid$ideo> -i & ddid$ideo < i
  the.center <- subset(ddid,center==T)
  center <- as.data.frame(as.matrix(by(the.center$ss,the.center$year,sum)))
  names(center) <- paste("center",i,sep="_")
  all.centers <- merge(all.centers,center,by.x="year",by.y=0,all=T)
} 
all.centers$legis <- 48:57

pdf(file="fig_E1centerdefinitions.pdf",
    width=8,height=6)
par(mar=c(6,4,2,4))
plot(all.centers$year,all.centers$center_0.45,type="n",bty="n",xaxt="n",
     ylab=c("Centrist Share of Congress"),ylim=c(0,.8),xlab="")
axis(side=1,at=all.centers$year,labels=paste("Legis\n",all.centers$legis),padj=.5)
lines(all.centers$year,all.centers$center_0.2,col=gray(0.3),lty=2)
lines(all.centers$year,all.centers$center_0.25,lty=1)
lines(all.centers$year,all.centers$center_0.3,col=gray(0.4),lty=3)
lines(all.centers$year,all.centers$center_0.35,col=gray(0.5),lty=4)
lines(all.centers$year,all.centers$center_0.4,col=gray(0.7),lty=1)
legend(x="topleft",bty="n",
       legend=c("[-.20,.20]",
                "[-.25,.25]",
                "[-.30,.30]",
                "[-.35,.35]",
                "[-.40,.40]"),
       lty=c(2,1,3,4,1),
       col=c(gray(.3),1,gray(.4),gray(.5),gray(.7)))
dev.off()



## Figure D1
## Reads PREPPS vs 2 data and compares it with our own estimates
## PREPPS conducted in 2018 and 2019, in between BLS9 and BLS9
## PREPPS data are available from http://wiesehomeier.net/prepps/

## 'PREPPS Latam V2.dta' was collected in 2018/2019 -> BLS8 (2017) and 9 (2021)
## 2015_LA_sum.dta was collected in 2015 -> BLS 8 (2017)
## Pilot_LA_sum.dta was collected in 2011 -> BLS 7 (2012)
## 2006_LA_sum.dta was collected in 2006/07 -> BLS 5 (2005) and 6 (2009)


## Get BLS estimates 
load("bls9_estimates_partiespresidents_long.RData") 
bls <- long.table
rm(long.table)
bls5 <- subset(bls,year==2005)
bls6 <- subset(bls,year==2009)
bls7 <- subset(bls,year==2013)
bls8 <- subset(bls,year==2017)
bls9 <- subset(bls,year==2021)

## Get PREPPS 2018/19  and merge with BLS8/9
library(haven)
prepps.2018 <- read_dta("PREPPS Latam V2.dta")
prepps.2018 <- subset(prepps.2018,dimensionname=="Left-Right"&
                   country=="Brazil")
prepps.2018$party.or.pres <- ifelse(prepps.2018$partyname_es=="",
                                    prepps.2018$presidentname,
                                    prepps.2018$partyname_es)
prepps.2018$party.or.pres <- car::recode(prepps.2018$party.or.pres,
                                    "'Partido Democrático Trabalhista'='PDT';
                                     'Democrátas'='DEM';
                                     'Partido da República'='PL';
                                     'Movimento Democrático Brasileiro'='MDB';
                                     'Progressistas'='PP';
                                     'Partido Socialista Brasileiro'='PSB';
                                     'Partido da Social Democracia Brasileira'='PSDB';
                                     'Partido dos Trabalhadores'='PT';
                                     'Partido Trabalhista Brasileiro'='PTB';
                                     'Partido Republicano Brasiliero'='REP';
                                     'Partido Social Democrático'='PSD';
                                     'Podemos'='PODE';
                                     'Partido Novo'='NOVO';
                                     'Partido Socialismo e Liberdade'='PSOL';
                                     'Partido Verde'='PV';
                                     'Rede Sustentabilidade'='REDE';
                                     'Partido Social Liberal'='PSL';
                                     'Michel Temer'='TEMER';
                                     'Jair Bolsonaro'='BOLSONARO'")

## Merge estimates ##
d <- merge(subset(bls8,select=c(party.or.pres,ideo,ideo.se)),
           subset(bls9,select=c(party.or.pres,ideo,ideo.se)),
           by="party.or.pres",all=T,suffixes=c("8","9"))
# Average BLS8 and 9 estimates to simplify comparisons
d$ideo <- apply(subset(d,select=c(ideo8,ideo9)),1,mean,na.rm=T)

d.2018 <- merge(subset(d,select=c(party.or.pres,ideo)),
           subset(prepps.2018,select=c(party.or.pres,score_pos)),
           by="party.or.pres",all=T)
d.2018$year <- 2018

# Get PREPPS 2015 and merge with BLS8
# General left-right dimension is dimension 13 in the dataset
prepps.2015 <- read_dta("2015_LA_sum.dta")
prepps.2015 <- subset(prepps.2015,Dimension==13&Country=="BRA"&
                        Party!="Barack Obama"&Party!="Evo Morales"&Party!="Juan Manuel Santos")
prepps.2015 <- prepps.2015 %>%
  select(Party,Mean) %>%
  rename(party.or.pres=Party,score_pos=Mean) %>%
  mutate(party.or.pres=gsub(" ","",toupper(party.or.pres)))
prepps.2015$party.or.pres <- car::recode(prepps.2015$party.or.pres,
                                         "'PMDB'='MDB';
                                     'PPS'='CID';
                                     'PR'='PL';
                                     'PRB'='REP'")  
## Merge estimates ##
d.2015 <- merge(subset(bls8,select=c(party.or.pres,ideo))
                ,prepps.2015,by="party.or.pres",all=T)
d.2015$year <- 2015

# Get PREPPS 2011 and merge with BLS6
# General left-right dimension is dimension 13 in the dataset
prepps.2011 <- read_dta("Pilot_LA_sum.dta")
prepps.2011 <- subset(prepps.2011,Dimension==13&Country=="Brazil")
prepps.2011 <- prepps.2011 %>%
  select(PartyName,Mean) %>%
  rename(party.or.pres=PartyName,score_pos=Mean) 
prepps.2011$party.or.pres <- car::recode(prepps.2011$party.or.pres,
                                         "'Partido Democrático Trabalhista'='PDT';            
                                          'Democrátas'='DEM';                                 
                                          'Partido da República'='PL';                      
                                          'Partido do Movimento Democrático Brasileiro'='MDB';
                                          'Partido Progressista'='PP';                      
                                          'Partido Socialista Brasileiro'='PSB';              
                                          'Partido da Social Democracia Brasileira'='PSDB';    
                                          'Partido dos Trabalhadores'='PT';                  
                                          'Partido Trabalhista Brasileiro'='PTB';             
                                          'Dilma Rousseff'='DILMA';                             
                                          'Partido Verde'='PV';                             
                                          'Partido Popular Socialista'='CID';                 
                                          'Partido Social Cristão'='PSC';
                                          'Partido Socialismo E Liberdade'='PSOL';
                                          'Partido Comunista do Brasil'='PCDOB'")  

## Merge estimates ##
d.2011 <- merge(subset(bls7,select=c(party.or.pres,ideo))
                ,prepps.2015,by="party.or.pres",all=T)
d.2011$year <- 2011

# Get PREPPS 2006/07 and merge with BLS5 and 6
# General left-right dimension is dimension 13 in the dataset
prepps.2006<- read_dta("2006_LA_sum.dta")
prepps.2006 <- subset(prepps.2006,Dimension==13&Country=="BRA")
prepps.2006 <- prepps.2006 %>%
                select(Party,Mean) %>%
                rename(party.or.pres=Party,score_pos=Mean) %>%
                mutate(party.or.pres=gsub(" ","",toupper(party.or.pres)))
prepps.2006$party.or.pres <- car::recode(prepps.2006$party.or.pres,
                                         "'PRES'='LULA';
                                          'PPS'='CID';
                                          'PFL'='DEM';
                                          'PMDB'='MDB'")  

## Merge estimates ##
lula <- bls[min(which(bls$party.or.pres=="LULA")),]
bls5 <- rbind(bls5,lula) #add Lula's retrospective rating here
d <- merge(subset(bls5,select=c(party.or.pres,ideo,ideo.se)),
           subset(bls6,select=c(party.or.pres,ideo,ideo.se)),
           by="party.or.pres",all=T,suffixes=c("5","6"))

# Average BLS5 and 6 estimates to simplify comparisons
d$ideo <- apply(subset(d,select=c(ideo5,ideo6)),1,mean,na.rm=T)
# Merge BLS and PREPPS
d.2006 <- merge(subset(d,select=c(party.or.pres,ideo)),
                subset(prepps.2006,select=c(party.or.pres,score_pos)),
                by="party.or.pres",all=T)
d.2006$year <- 2006



### Stack all merged waves
d <- na.omit(rbind(d.2006,d.2011,d.2015,d.2018))


# Correlation reported in main body of the research note:
cor.test(d$ideo,d$score_pos)

# Figure in appendix  
pdf(file="fig-D1blsprepps.pdf",
    width=8,height=8)
par(mar=c(4,4,1,4))
pres <- ifelse(is.element(d$party.or.pres,c("LULA","DILMA","TEMER","BOLSONARO")),T,F)
plot(d$ideo,d$score_pos,type="n",bty="n", 
     xlab="BLS Position Estimates",
     ylab="PREPPS Overall Left-Right Score",xlim=c(-1.2,1.2)) 
mtext(expression("" %<-% Left),side=1,line=1, adj = 0,cex=1)
mtext(expression(Right %->% ""),side=1,line=1, adj = 1,cex=1)
mtext(expression("" %<-% Left),side=2,line=1, adj = 0,cex=1)
mtext(expression(Right %->% ""),side=2,line=1, adj = 1,cex=1)
the.reg <- lm(d$score_pos~d$ideo)
abline(reg=the.reg,col=gray(0.6))
points(d$ideo,d$score_pos,
       pch=ifelse(pres,24,21),bg=ifelse(pres,"white",gray(0.6)))
legend(x="bottomright",
       legend=c("Party","   PT","   PSDB","President"),
       pch=c(21,21,21,24),
       pt.bg=c(gray(0.6),gray(0.6),gray(0.6),"white"),
       col=c(1,1),
       bty="n")
legend(x="bottomright",
       legend=c("Party","   PT","   PSDB","President"),
       pch=c(21,3,3,24),
       pt.bg=c(gray(0.6),1,"white","white"),
       col=c(1,1,"white",1),
       bty="n")
text(d[pres,"ideo"],d[pres,"score_pos"],
     labels=paste(d$party.or.pres[pres],c("","I","II","","")),
     pos=c(4,4,1,4,1),xpd=NA,cex=0.8)
rede <- which(d$party.or.pres=="REDE")
text(d$ideo[rede],d$score_pos[rede],labels="REDE'18",pos=3,cex=0.9)
mdb <- which(d$party.or.pres=="MDB"&d$year==2015)
text(d$ideo[mdb],d$score_pos[mdb],labels="MDB'15",pos=4,cex=0.9)
#dem <- which(d$party.or.pres=="DEM")
#text(d$ideo[dem],d$score_pos[dem],labels="DEM",pos=3,cex=0.9)
pt <- which(d$party.or.pres=="PT")
points(d$ideo[pt],d$score_pos[pt],pch=3)
#text(d$ideo[pt],d$score_pos[pt],labels="PT",pos=3,cex=0.9)
psdb <- which(d$party.or.pres=="PSDB")
points(d$ideo[psdb],d$score_pos[psdb],pch=3,col="white")
dev.off()


